Dynamical Modularity in Automata Models of Biochemical Networks

Given the large size and complexity of most biochemical regulation and signaling networks, there is a non-trivial relationship between the micro-level logic of component interactions and the observed macro-dynamics. Here we address this issue by formalizing the concept of pathway modules developed by Marques-Pita and Rocha [1], which are sequences of state updates that are guaranteed to occur (barring outside interference) in the causal dynamics of automata networks after the perturbation of a subset of driver nodes. We present a novel algorithm to automatically extract pathway modules from networks and characterize the interactions that may take place between the modules. This methodology uses only the causal logic of individual node variables (micro-dynamics) without the need to compute the dynamical landscape of the networks (macro-dynamics). Specifically, we identify complex modules, which maximize pathway length and require synergy between their components. This allows us to propose a new take on dynamical modularity that partitions complex networks into causal pathways of variables that are guaranteed to transition to specific dynamical states given a perturbation to a set of driver nodes. Thus, the same node variable can take part in distinct modules depending on the state it takes. Our measure of dynamical modularity of a network is then inversely proportional to the overlap among complex modules and maximal when complex modules are completely decouplable from one another in the network dynamics. We estimate dynamical modularity for several genetic regulatory networks, including the full Drosophila melanogaster segment-polarity network. We discuss how identifying complex modules and the dynamical modularity portrait of networks explains the macro-dynamics of biological networks, such as uncovering the (more or less) decouplable building blocks of emergent computation (or collective behavior) in biochemical regulation and signaling.

1 introduction 3 individual node interactions give rise to global properties of the network such as its modularity or attractor landscape. It is, therefore, useful to understand how the micro-dynamics of component interactions give rise to the macro-dynamic behavior observed in networks.
Modularity methods have become a standard approach to decompose graphs and solve problems such as community detection. These approaches typically focus on topological modules, where a module is defined by a higher density of links between nodes in the same module and a lower density of links to nodes in different modules [17,18]. Such structural modularity measures ignore the dynamical information of the network and may, therefore, be limited when it comes to discovering functional modules [19,20]. It has also been shown that relying on structural information alone cannot account well for dynamics in control problems on BNs [21].
Additionally, the same circuit topologies in biological systems can switch between distinct steady states and dynamical behaviors (multi-stability and multifunctionality), which indicates that there is not a clear relationship between structure and function [22]. show that controlling these motifs can drive a network to a desired attractor [27].
We use the threshold network representation of automata network dynamics proposed by Marques-Pita and Rocha, the dynamics canalization map (DCM), which represents both Boolean states of every node. They simplified the description of a network's dynamics by explicitly removing redundancy in the logical expressions governing each node's update (described in the following subsections). This redescription allows to describe dynamical modularity, critical nodes that guarantee convergence to steady states (with incomplete knowledge of initial conditions), and measures of macro-level canalization [1]. Using the DCM, they identified the dynamical modules on the Drosophila melanogaster single-cell segment polarity network (SPN) [28] that are controlled by the network's inputs.
Here we formally define modules on the DCM. We expand the analysis of [1] by finding such modules not just for the inputs but for all possible nodes and low-order interactions on the drosophila segment polarity network (SPN). We develop theoretical results to explain the interaction of different dynamical modules and define synergy between modules. The concept of synergy is used to define a subset of modules, called complex modules, that capture all unique dynamical information in the BN.
We then use complex modules to define a measure of dynamical modularity for BNs. This measure is driven purely by the network's dynamics rather than its structure. Unlike other measures, our dynamical modularity measure differentiates between different states of a variable. It accounts for the fact that a single variable may be part of multiple functional modules depending on its behavior and thus may be useful in models of multifunctionality in biological systems. Additionally, our measure makes explicit the influence of each variable in each state and the control that a given variable has on the rest of the network. We also use the mean dynamical modularity of a network to quantify its decomposability and use this metric to compare different biological networks. This paper is organized as follows: In the remainder of the Introduction, we give some background on schematic redescription, canalyzing maps, and the DCM introduced in [1]. In section II, we give theoretical results concerning the formal definition of pathway modules (based on the perturbation of seed nodes), complex modules (which have synergy between seed nodes and are of maximal size), and dynamical modularity (based on an optimal partition of the DCM). In section III, we apply our methodology to the drosophila single-cell and parasegment SPN to find complex modules on these networks, quantify their dynamical modularity, and compare them to other biological networks of comparable size. In section IV, we offer our conclusions and note some limitations to our methodology.

Boolean Networks
A Boolean automaton x can be in one of two states, x ∈ {0, 1} (traditionally ON or OFF, representing activation of the molecule of interest above or below a certain threshold) [14,15,16]. The state of an automaton is updated in discrete time steps governed by a logical function of its inputs (e.g., A = B ∨ C). More formally, x t+1 = f(x t 1 ,...,x t k ) for automaton x with k inputs at time step t, where the mapping f : {0, 1} k → {0, 1} can be read from the automaton's look-up table F that denotes the output x t+1 for each of its 2 k input combinations (refer to Figure 1).
A BN is a graph B = (X, E) where X is a set of n Boolean automata nodes and E is a set of directed edges e i,j ∈ E that indicate that node x i is an input to node x j [14,29]. The configuration of a BN is the value of all automata states x ∈ X. Automata nodes within the network may be updated synchronously in discrete time steps or asynchronously, where each node has a different update schedule (deterministic asynchronous BNs) or is selected at random with some probability [30,29]. Deterministic BNs are guaranteed to eventually resolve in some attractor, either a fixed point or a limit cycle of repeating states. In biological networks, attractors are typically associated with some phenotype (such as a wildtype expression pattern) [16].

Node canalization
Canalization refers to the ability of a subset of inputs to control an automaton's transition [31,32,33,34,1,21,35]. The Boolean transition function of the automaton is thus satisfied by this subset of inputs and the logical states of its other inputs can be ignored. This shows inherent redundancy in the transition function that can be removed via schematic redescription, a process introduced in [1]. The Quine-McCluskey Boolean minimization algorithm is used to reduce an automaton's look-up table F to a set of prime implicants [36] that are then combined to form a set of wildcard schemata F .
Within the wildcard schemata only essential inputs (enputs) are defined as ON or OFF. Every other input is replaced with the wildcard symbol (#) which indicates that those inputs are redundant (given the states of the enputs) and do not affect the automaton's state transition [1]. Additionally, permutations of inputs that leave the transition unchanged are marked with a position-free symbol ( o ) and combined to form group-invariant enputs. This captures the input-symmetry within the wildcard schemata and further reduces them to a set of two-symbol schemata without permutation redundancy F (see Figure   1). It is not surprising to find input redundancy and input-symmetry at the automaton level (known as micro-canalization) in biological networks because they are known to be robust, which buffers perturbations to individual inputs and can help maintain wildtype phenotypes [37,33,10].
After all redundancy is removed, the necessary and sufficient logic of automata nodes can also be represented by a canalyzing map (CM) [1], a type of threshold network [38]. The CM of a variable x represents all logical transitions (two-symbol schemata) that can result in an update of the state of x. It is composed of s-units which represent nodes (x and its inputs) in a particular discrete state (e.g., ON or OFF), t-units that implement the transition function of x via thresholds, and fibres that connect s-units and t-units together. A fibre may connect one s-unit to one t-unit (indicating that the s-unit sends one signal to the t-unit), or multiple fibres may be merged together so that several s-units send one signal to a t-unit (representing permutation redundancy). T-units have a threshold value τ and only fire if they receive at least τ simultaneous incoming signals from s-units; s-units fire if they receive any incoming signal from a t-unit. All logical interactions between x and its inputs are represented such that each t-unit corresponds to one schema in F . Thus the CM is equivalent to the schemata necessary to ensure the logical transition from x t to x t+1 (see Figure 1). The node x has two inputs. Center: The look-up table F describing the transition of x and the one and two-symbol schemata redescription (F and F ; # is a wildcard enput, o is a group-invariant enput). Right: The CM for x with the associated logical function representing both state updates. Black s-units represent nodes in their ON state, white s-units represent nodes in their OFF state, and blue diamonds represent t-units labeled with their threshold value (τ). Each edge (fibre) represents the logical contribution of a source unit to an end unit and thus dictates how X will turn ON or OFF. The CM shows that one ON input is sufficient to turn X ON but that two OFF inputs are needed to turn X OFF.

Dynamics canalization maps
The canalyzing map of each automaton can be be integrated into a larger threshold network called the DCM [1]. This parsimonious network, with input and symmetry redundancy removed, represents the control logic and dynamics of the entire BN. Each possible state of each variable is represented in the DCM; thus there are 2n s-units for a BN of size n and t-units that represent each schema, as defined for each individual CM.
Because the DCM captures the state updates of all automata in the network, it is possible to determine a logical sequence of signals based on the deterministic firing of s-units and t-units. Given an initial set of s-units (and barring any outside interference), a logical sequence of state updates is guaranteed to occur based on the logic embedded in the schema F and represented by the t-units in the DCM. This sequence of state updates is referred to as a pathway module because it is a sequence of s-units and t-units firing in the DCM and the state updates occur independently of all other node states that are not involved in the firing of t-units [1]. Importantly, the initial conditions of a pathway module may only define the state of a (small) subset of nodes in the network, while the complement of this set is assumed to be in an unknown state (represented by '#'). These known and unknown node states are referred to as a partial configurationx of the network. Because dynamics are deterministic, this partial configuration causes a sequence of state updates (which we refer to as dynamical unfolding) that results in an outcome configuration P, which may be either another partial configuration or an attractor of the network. If P is an attractor then we know thatx controls the full network dynamics because the attractor is guaranteed to occur based on the subset of known s-units in the initial condition if the initial node states are sustained [1]. Thus, we can compute global (macro-scale) dynamical information from only partial knowledge of the network configuration by integrating knowledge of the local (micro-scale) dynamics.
We note that the calculation of dynamical unfolding is similar to the calculation of partial steady states in the logical interaction hypergraph [39], cascading effects of node removal and the calculation of elementary signaling modes on the expanded network [40], and calculation of the logical domain of influence on the expanded network [41]. However, none of these methods explicitly remove input symmetry from the transition functions of nodes or consider different types of node perturbation.
In the next section we give a formal treatment of pathway modules and provide definitions of novel concepts such as pathway module interaction, complex modules, and measures of dynamical modularity.

formal results
Pathway modules can be discovered using a breadth-first search algorithm on the DCM starting from the initial seeds. However, unlike a search on a traditional graph, the different types of units on the DCM must be considered.
For a t-unit to be included, its threshold τ must be met. Furthermore, the threshold is met only if its inputs fire simultaneously; therefore, the firing of a t-unit with multiple inputs requires confidence that all of its inputs can fire at the same time. If the initial seeds are pinned (not allowed to change), then they fire at every time step and all downstream s-units fire repeatedly as well.
However, if the initial seeds are only perturbed once (pulse perturbation), (a) (b) Figure 2: Example GRN. The interaction graph is shown on the left, the DCM on the right. This network has two inputs (i1,i2), two genes (g1,g2), and two proteins (P1,P2 that modules based on pinning perturbation may uncover fixed points in the network but not limit cycles.
In the example GRN in Fig Multiple pathway modules may simultaneously unfold within the complex dynamics of a network. Thus, it is natural to explore how different pathway modules may interact, assuming that their initial conditions are simultaneously activated. This leads to concepts such as independence, logical obstruction, and synergy. Synergy is used to define complex modules, which are maximal pathway modules whose seeds have synergy between them; complex modules are further used to define the dynamical modularity of the network. Formal definitions of these concepts are given in the following subsections.

Defining Pathway Modules
The definition of the DCM is provided in the Introduction. There is an important constraint on the set S t of s-units that fire at time t. If an s-unit that represents the state of x fires at time t, the s-unit that represents the negation of the state of x cannot fire. Thus, x-1 ∈ S t =⇒ x-0 / ∈ S t and vice versa. S t i represents a partial configuration whereby the logical state of the corresponding node variables x ∈ X is known, but any other variable may be in either state (ON or OFF). The number of (full) configurations that a partial configuration S t describes is 2 |X|−|S t | . As a corollary, max(|S t i |) = |X| denotes the situation when the logical state of all variables is known, and we know the precise configuration the network is in A seed set is a set of s-units S 0 ⊂ S that, without loss of generality, are assumed to fire at time t = 0. This functions as an external perturbation applied to a subset of variables X 0 ⊂ X of network B. Note that these s-units can represent variables in either the ON or OFF states. S 0 can be treated differently in the dynamical unfolding process. When studying pulse perturbations, s-units s ∈ S 0 are assumed to fire only at t = 0, unless the logic of the canalized dynamics (µ) leads them to fire again. In contrast, pinning perturbations lock the s-units s ∈ S 0 into firing at every time step t 2 . In this case, we always have S t ⊂ S t+1 in the canalized dynamics.
Unless otherwise noted, we use pinning perturbation, which represents situations where the perturbation signal is assumed to be in steady-state and last much longer than the transient dynamics that it affects. This is a realistic scenario, especially for gene regulatory network models.
A module in the DCM is defined as M t ≡ (S t , θ t ). It is a tuple consisting of a set of s-units S t ⊂ S and t-units θ t ⊂ Θ that fire at time t due to a perturbation applied to seed set S 0 at t = 0. is T , whereas its size is simply the cardinality of its corresponding pathway module set |S|. Finally, P s denotes the set of all pathway modules M(S 0 ) that unfold in a DCM from seed sets of size s = |S 0 |.

Interaction of pathway modules
Given that there is an exponential number of pathway modules, for computational purposes we must narrow down the possible dynamical interactions of the system. Therefore, we define how pathway modules may interact to hone in on those that best define the computational dynamics of the network. Specifically, we consider that two pathway modules M i and M j can be combined by considering the dynamical unfolding of the union of their seed sets, for all t T (the sequence is 3 We assume that external signals take precedence over natural dynamics of the network. As such, a perturbation that causes an s-unit to fire will prevent any logically contradictory s-units from firing during the same time step. For pinning perturbation, a variable may only fire in one state during the dynamical unfolding of a pathway module; that is, x-1 ∈ S =⇒ x-0 / ∈ S and vice versa. Thus, pinning perturbation cannot lead to limit cycles or other oscillations. 4 Seed sets can, of course, be combined via set-theoretical operations at distinct time steps, but here we only consider unions of seed sets at time t = 0. preserved temporally, where constant k 0 indicates that the dynamical unfolding of the subsequence M i within M j begins k time steps after t = 0).
In cases where there is also no synergy between M i and M j (see below), . Logical obstruction can only occur when the union of the seed sets results in a logical contradiction whereby s-units for the state of at least one node variable and its negation are both included (in the seed sets or their dynamical unfolding logic), for example, x-1 ∈ S 0 i ∧ x-0 ∈ S 0 j . This means that combining the pathway modules results in one or more s-units not firing that would have fired had the pathway modules been separate.
Synergy exists between pathway modules M i and M j if their combination M i,j results in an s-unit (or t-unit) firing during time step t = k when that s-unit (or t-unit) does not fire in either pathway module M i or M j at time If there is no logical obstruction (see above) between M i and M j , then: . This means that additional s-units fire as a result of combining the pathway modules.
Pathway modules M i and M j are decoupled if there is no partial subsumption (s-unit overlap), logical obstruction, or synergy between them. This | are necessary and sufficient conditions for pinning perturbation.
this means that all modules M i ∈ Λ 1 are complex. We define I s ⊂ P s as the set of all complex modules with seed size s.
To be a complex module then, a pathway module must be maximal and every submodule (seed combination) within it must add synergy. For a given s, each complex module M i ∈ I s represents the introduction of a unique firing sequence that is not part of any other pathway module with seed set size x s (otherwise it would not be maximal) or any combination of pathway modules with seed set sizes x s (because every submodule adds synergy). As a consequence, complex modules are the building blocks of a network's dynamics. Importantly, complex modules allow the vast number of possible dynamical interactions to be reduced to a smaller set of unique signaling sequences that may represent (full or partial) biologically functional pathways.
Given a DCM, its complex modules are defined using only parameter s, the size of the seed set. Moreover, the complex modules are hierarchical in the sense that modules with larger seed set sizes are composed of (partially or fully subsumed) modules with smaller seed set sizes. Higher-level complex modules thus capture the emergent behavior (e.g., synergy) that results from combining lower-level modules. For s = 1, the set of complex modules I 1 = Λ 1 is the minimal set of pathway modules that covers the DCM.
However, synergy allows for the DCM to be covered using fewer pathway modules (with larger seed sets).

Dynamical Modularity
With the characterization of pathway module interaction and the introduction of complex modules, we may now characterize the macro-scale properties of the network, such as its organization and modularity. The concept of dynamical modularity is analogous to structural modularity in graphs, whereby it is considered to occur when interactions within modules are strong and interactions between modules are weak. Dynamical modularity is, however, a more general concept derived from Simon's framing of complex systems in terms of near-decomposability [45], which includes both structural and dynamical interactions [46, 24, 21]. We use it to find subsets of variables and their dynamical states that are more or less easy to decouple in a given multivariate dynamical system. Thus, we can compare complex (dynamical) networks (or systems) to determine which ones are easier to decouple into independently-functioning parts.
The cover of a DCM is a set of pathway modules Π = {M i , M j , ...} such that the union of their corresponding pathway modules sets is equivalent to the set of all s-units in the DCM; that is, S Π = i S i = S for set S i with corresponding pathway module M i ∈ Π. Because pathway module sets may overlap, this results in a fuzzy partition of s-units. |Π| is the size of the cover in terms of the number of pathway modules that comprise it. In practice, we are interested in finding a cover Π s that is composed of modules with no more than s seeds; additionally, we focus on small s because modules with smaller seed sets are generally easier to control. The independence of a pathway module M i from a set of pathway modules Σ is the fraction of unique s-units that fire within the pathway module: gives us the independence of pathway module M i within a cover Π.
The dynamical modularity of a cover Π is the distribution of independence scores for all its pathway modules: can derive several meaningful statistics from this distribution, but here we use its mean value as a characterization of dynamical modularity: As there are many possible covers of the DCM, it is useful to find one that is optimal for measuring modularity. One can define optimality of a cover in different ways, for example, the cover with the minimal size, Π min . However, as the goal is to find modules that are decoupable from one another, here we define optimality by maximizing the mean dynamical modularity D(Π).
We define Π * s as the optimal cover of the DCM, where Π * s has the maximum mean dynamical modularity among covers that satisfy the property that M i ∈ I for all pathway modules M i ∈ Π s (i.e., the cover is composed only of complex modules with seed set size |S 0 | s). The corresponding mean dynamical modularity of this cover indicates how decouplable a network is into its constituent dynamical building blocks.
We note that the optimal cover may be difficult to calculate exactly due to the exponential number of possible module combinations. When networks and seed set sizes are too large to find an exact solution, we use a greedy selection to estimate Π * s . The algorithm starts with an empty set Σ = {} and then iteratively adds the complex module M i ∈ I s that has the highest independence score from Σ, max(ind(M i , Σ)), until a cover is reached, such that S Σ = S. We use this cover Σ to estimate the optimal cover and its mean dynamical modularity D(Σ) to estimate the optimal mean dynamical modularity of the network D(Π * s ). This quasi-optimal solution provides a lower-bound on D(Π * s ). In practice, multiple covers may have the same mean dynamical modularity, meaning that there may be multiple optimal solutions.
The modules that compose a cover Π s have a distribution of seed set sizes that are bounded by s. By definition, the mean dynamical modularity of Π * s should increase or remain the same as s increases, because each increase in maximum seed set size allows for more modules to choose from when finding an optimal cover. The characteristic seed number s * occurs when increasing s no longer increases the mean dynamical modularity of Π * s . This represents the minimal seed set size that allows for a network to optimally be dynamically decoupled.
In Fig. 4 panel f, the six complex modules shown comprise Λ 1 and the as- The mean dynamical modularity of this cover is D(Π * 1 ) = 0.71. Additionally, there are two complex modules (assuming pinning pertur- . With pinning perturbation, this network can be covered by three modules: with mean dynamical modularity d = 0.6. Note that this minimal cover is not optimal because D(Π * 1 ) has a higher mean dynamical modularity. Rather, the optimal cover at s = 2 is Π * The dynamical modularity does not increase with greater s so the characteristic seed number of this network is In summary, dynamical modularity allows for the comparison of dynamical organization between networks. Dynamically simple networks will be easier to decouple and will tend to have higher mean dynamical modularity and lower characteristic seed number. More complex networks will have more overlap dynamically and will tend to have lower mean dynamical modularity and higher characteristic seed number.
We note that, because our dynamical cover partitions s-units in the DCM, the same node will belong to different modules depending on what state it is in. This is a more granular but lengthier description of the dynamics of the system than considering the nodes themselves.
Computationally, our optimal measure of dynamical modularity is NPhard as it entails the set cover optimization problem [47,48] (by comparing all possible sets of complex modules that cover the set S). However, the number of modules to consider is greatly reduced by considering only complex modules rather than all pathway modules. For all real-world networks that we observed, I s << P s because many pathway modules are either subsumed into longer modules or they have non-synergistic seeds. Thus, by only considering complex modules, we can greatly reduce the computational complexity of our dynamical modularity optimization problem.

experimental results
We demonstrate our approach by using the above formalism to analyze the dynamical modularity of experimentally-validated models of biochemical regulation from systems biology. We focus on the single-cell and full as the single-cell model, plus one SLP input per cell (which represents the output of upstream developmental pathways and is usually assumed to be pinned ON or OFF). Furthermore, each of the 4 cells in the parasegment regulates its two neighboring cells via its internal WG and HH nodes, which act as external, inter-cellular inputs to the two neighboring cells. Thus, in the parasegment SPN there are no input nWG or nHH nodes for each cell, which results in a network of n = 60 variables. This model has been extensively studied [53,28,54,55,56,57,1,21] and its attractors are fully known for the single-cell case, which makes it a useful example to test our modularity methodology.
We compute all complex modules present in the dynamics of the singlecell and full parasegment SPN models for seed set size s 6. This analysis allows us to calculate the dynamical modularity of the SPN and compare it to the other systems biology models we analyzed. In the present work we only consider the dynamics that ensues from pinning perturbation of the seed set.
This form of perturbation is feasible for experimental work with the target biological models (e.g., gene silencing in drosophila [58]), but our approach also allows for other forms of perturbation (e.g., pulse perturbations) that can be explored in the future.

Complex modules in the Drosophila Single-Cell SPN
Analysis of the drosophila single-cell SPN highlights many useful features of our dynamical modularity methodology. It allows us to reduce the enormous complexity of this dynamical system's state space to only a few building blocks, characterize which dynamics are involved in transient configurations and attractors, understand the dynamical overlap between similar modules, and study the specific effects of individual nodes and sets of nodes (e.g., comparing network inputs to internal nodes).
To understand the effect of perturbing a specific node in the single-cell SPN with a brute-force approach would require analysis of 2 17 configurations.
Additionally, there are 2 s 17 s different ways to perturb s seeds, where each node x in the seed set of size s may be in the state x-0 or x-1. However, our dynamical modularity methodology offers a much more direct approach without enumerating all network configurations and by reducing the total possible number of seed sets to only those that are complex modules. This results in a total of 56 complex modules. However, this number can be further reduced using the maximal seed heuristic where modules are discarded if they contain seeds that do not initiate maximal pathway modules.
That is, we consider core complex modules to be those where each seed's pathway module is itself maximal, M x ∈ Λ 1 , ∀x ∈ S 0 . This reduction removes modules that would have been subsumed if not for contradiction. For that only one of these modules will ever be active at a time in the dynamical system. Because of this independence, these modules and modules of greater seed set size that subsume them appear in the optimal covers of the DCM for seed set sizes s = [1,4].
Finally, these complex modules help to elucidate the role of specific nodes, such as the inputs, in the single-cell SPN dynamics. It is already known that the input nodes control much of the network dynamics [1,21]; however, our analysis shows when the inputs themselves may result in an attractor and when an additional node such as wg or PT C is needed 6 . Additionally, we see the effects that internal nodes, such as PT C or CIR, have by themselves. This is especially useful in analyzing the parasegment SPN (n = 60) because all nodes in this network are internal except for the four SLP inputs. maintains the state of its inputs (and, therefore, its downstream products) due to a positive feedback loop between PT C-1 and nHH-0.

Complex modules in the Drosophila Parasegment SPN
We use our dynamical modularity methodology to aid the analysis of the drosophila parasegment SPN. As in the single-cell case, this allows for a great reduction in the complexity of the dynamical system model and enables us to characterize dynamical trajectories within the system. Using complex modules, we are able to find a minimal seed set that reaches the wildtype attractor of drosophila that is lower than previous estimates seen in the literature. Furthermore, we are able to characterize the transient dynamics toward this attractor in terms of the intracellular modules we found in the single-cell SPN.
It is computationally restrictive to discover all complex modules for the drosophila parasegment model, but as in the single-cell case, the number of complex modules found for low seed set sizes is substantially less than the number of possible pinned perturbations. We again use the maximal seed heuristic to speed up computation time. For s = 1, we find 40 complex modules out of 120 possible pathway modules; for s = 2, we find 40 core complex modules out of 7080 pathway modules; for s = 3, we find 152 core complex modules out of 273760; for s = 6, we find 1475 out of more than 3.2 billion. Despite this great reduction, there are still enough core complex modules to make analysis difficult; therefore, we focus only on the modules with the greatest size for each s value. We are aided also by the symmetry in the model (the index of each of the four cells does not matter because each one is identical to the next) and the fact that the intracellular dynamics of each cell in the parasegment is the exact same as in the single-cell case, except that they are influenced by inputs from neighboring cells (WG and HH) as well as the external input SLP. for s = 6, the largest complex modules resolve all 60 node states, which indicates a network attractor (see supplementary materials). By sampling larger pathway modules, we find that 9 nodes are sufficient to reach the wildtype attractor of the parasegment SPN, which is lower than previous estimates seen in the literature [1,44,59] (see Fig. 6).
The discovery of the complex modules allows us to describe the transient dynamics that occur in the parasegment model toward the wildtype attractor (see Fig. 6 and supplementary materials) and other attractors in the network.
We find that inter-cellular interactions can be described by the intracellular complex modules found for the single-cell SPN, where the unfolding of one intracellular module M i may initiate another module M j by causing the seed of that module to fire (S 0 j ∈ S i ). We observe that these transient dynamics move back and forth across cell boundaries as the unfolding of one intracellular module affects its neighboring cells (see the example of the wildtype attractor in Fig. 6 and further examples in supplementary materials).
For example, a Σ 2 module will initiate M nHH−1 in both adjacent cells (which will interact synergistically with any Σ 1 or Σ 2 modules already active in those cells) and Σ 2 modules in non-adjacent cells will activate Σ 3 modules in their neighboring cells (via neighboring WG-0 and HH-1). Interestingly, some intracellular modules that were dynamically independent in the single-cell SPN work together across cell boundaries in the parasegment. For example, M (SLP−1,nHH−1) (a Σ 3 module) will activate WG-1, which can interact with SLP-0 in an adjacent cell to initiate M nWG−1,SLP−0 (a Σ 2 module) in that cell, even though Σ 2 and Σ 3 modules are dynamically independent (with the exception of CIR-0) within the same cell.
By analyzing the largest complex modules, we find that internal nodes play a large role in control of the dynamics, particularly en, wg, PT C, and We observe that the unfolding of pathway modules indicates a temporal order of events, such that one intracellular module may need to completely unfold to initiate another module in the same cell or an adjacent cell. Of course, during the natural dynamics of the network, multiple modules will be initiated simultaneously based on initial conditions and this natural activation may speed up convergence to a final state. However, the unfolding of a pathway module M i guarantees that any downstream node states x ∈ S i will occur (and gives an upper bound on the number of iterations it will take to do so via the module length), as long as the seed nodes are pinned for a long enough period. As an example, pinning the 9 seeds indicated in Fig.   6 (4 of which are the SLP inputs) is guaranteed to drive the network to the wildtype attractor in no more than 16 iterations. Given that this is a steady state, 16 iterations is also an upper bound on how long the seed nodes must be controlled, as the network will remain in the wildtype configuration once it reaches it.
As seen above, using the easily discoverable intracellular complex modules allows for a concise description and qualitative understanding of what is happening in the parasegment network, even though the network is too large for an exhaustive description. Because the drosophila parasegment is divided into four identical and modular cells, we can understand the dynamics of the larger network by understanding the complex modules within the individual cells and how those modules interact across cell boundaries. This gives us a clearer picture of the multivariate dynamics in the network than we could get via simple attractor analysis, while at the same time (due to its Boolean nature) being easier to analyze than a more detailed model (e.g., parameter estimation of differential equations). To better understand how signals propagate between dynamically dependent modules, we simplify the macro-level dynamics of the system and, thus, how the network computes.

Dynamical Modularity of the Drosophila SPN
We use the complex modules found in the drosophila SPN to estimate dynamical modularity of the single-cell and parasegment models by finding an optimal cover of the DCM, as described in section 2.3. For seed set size s = 1, the optimal cover is defined by the set of complex modules Λ 1 ; however, the general problem of finding an optimal cover composed of modules with at most s seeds is NP-hard based on the set-cover optimization problem [47,48]. We can, however, restrict the number of modules in the cover to a maximum defined by parameter q. For I complex modules there are, therefore, I q options to consider for the optimal cover. In general, dynamical modularity scores increase with q because there are a greater number of module combinations to choose from. However, for larger networks like the parasegment SPN, it is computationally prohibitive to estimate covers with even a moderately high q. In this case, we estimate the cover using a greedy algorithm (see section 2.3), where at each step the module with the highest independence score is added to the cover until q selections have been made.
In both cases, we only consider core complex modules using the maximal seed heuristic for inclusion in the cover.
For the single-cell SPN, the optimal cover at s = 1 is composed of 14 complex modules and has mean dynamical modularity D(Π * 1 ) = 0.69. For s = 2, we calculate the optimal cover Π * 2 for q 11; for s = 3, we calculate 1 2 3 4 Figure 6: Dynamical unfolding of the wildtype attractor in the drosophila parasegment. The minimal steady-state perturbation that we found sufficient to fully resolve the wildtype attractor is shown. Cell boundaries are separated by a blue line; white indicates that a node is OFF in that time step, black indicates that a node is ON, and grey indicates that the node state is unknown. The color of the node label indicates an associated complex module. Select modules are also highlighted in the dynamical unfolding process by colored borders. Arrows indicate the initiation of a new module; arrows pointing to the side or upward indicate that the respective s-unit is initiating a module in a neighboring cell. the optimal cover Π * 3 for q 8; for s 4 there are no more core complex modules so the solution does not change. We find that the optimal covers have more than the minimal number of modules necessary to cover the DCM for s > 1 (see Fig. 7, a-b). For s = 2, the DCM can be covered with 5 modules, but the optimal cover is |Π * 2 | = 10, D(Π * 2 ) = 0.73. For s = 3, the minimal cover is only three modules, but the optimal cover is |Π * 3 | = 8, D(Π * 3 ) = 0.81 (s * = 3 is also the characteristic seed number of the single-cell SPN). We also find that the DCM cannot be covered by only pinning the s-units representing the input nodes (SLP, nWG, nHH) because none of these modules include the s-units SMO-0 or CIR-1.
Comparing drosophila to other small GRNs, the yeast cell-cycle network We can also identify the modules that are included in the optimal cover.
For the parasegment SPN, the majority of the modules in our estimated covers are also of seed set size s = 1 with only a few larger complex modules.
Additionally, we see in Fig. 7d that the estimated mean dynamical modularity is roughly constant for both the parasegment SPN and the T-LGL leukemia network as s is increased. This suggests that even though the optimal cover is impossible to calculate at higher s values, and larger complex modules that include useful synergies may be missed, the mean dynamical modularity at s = 1 is a useful lower bound that may estimate well the true value for higher seed set sizes and may be used for comparison to other networks. The relative cover size |Π s |/2n is shown for each network's optimal cover Π * (based on maximizing D(Π s ), shown as full lines) and minimally-sized cover Π min (dashed lines) per maximum seed set size s. Note that |Π s |/2n 1 because there are 2n s-units that together create a trivial cover (every s-unit acts as a module's seed). The intermediate line |Π s | = n is indicated by grey dashes. We see that the optimal cover has more than the minimal number of modules for s > 1. (b) Same as panel a, but the dynamical modularity score d = D(Π s ) is shown for each network's optimal cover Π * (full lines) and minimally-sized cover Π min (dashed lines). (c) Results for the drosophila parasegment SPN and the T-LGL leukemia network. We estimate the optimal cover Π * s by greedy selection of the maximum module independence, max(ind(M i , Π s − M i )), as described in Section 2.3. The relative cover size |Π * s |/2n is shown for each network. Grey dashes indicate the line |Π s | = n. (d) Same as panel c, but the dynamical modularity score d = D(Π * s ) is shown for each network. Due to the suboptimal nature of the greedy selection, d is sometimes lower for larger maximum seed set sizes.

conclusion
We have formally defined the pathway modules first discovered in [1] and offered an algorithm to discover them via a modified breadth-first search on a network's DCM. We have also formally analyzed interactions between pathway modules and defined complex modules, which are maximal pathway modules with synergy (dynamic dependence) between seed s-units. This enables us to define the optimal cover of a DCM (a partition of s-units based on complex modules) and the associated mean dynamical modularity.
This modularity methodology is advantageous in that it greatly reduces the combinatoric number of pathway modules in Boolean GRNs for a given seed set size, 2 s n s , and only considers those that are complex modules. It also allows a bottom-up approach to defining modularity based on a single parameter, the maximum seed set size s. We were thus able to analyze all complex modules that occur in the drosophila single-cell SPN, and sample many that occur in the parasegment. This analysis of dynamical modularity allows an in-depth understanding of how dynamic signals propagate between cells in the parasegment and how dynamics unfold toward biologicallyrelevant attractors. It also provides a means to compare different GRNs based on the size and dynamical modularity of the optimal cover of their DCM for a given seed set size. This method can be applied to any automata network DCM, including those outside of the biological domain.
Importantly, this analysis describes the emergent computation that takes place in the example models as single seeds synergistically interact, causing their combined dynamical influence to logically propagate downstream, resulting in either an attractor or a partial configuration that represents the (much reduced) possible network configurations remaining. These partial configurations are useful for control problems because the description of pathway modules elucidates exactly which downstream nodes are guaranteed to be affected and which state they will be changed to. This can be especially helpful for target control [60, 41], where reaching a certain variable state is desirable (such as an apoptotic state in a cancer growth model) or not desirable (such as a proliferative state). Pathway modules, and complex modules in particular, can also be used to find seeds that have a similar effect on downstream targets. Alternative control strategies such as these have potential uses in biological networks, such as designing optimal therapeutic targets.
Pathway modules describe the network dynamics in a mechanistic, causal way. All s-units in a pathway module set are guaranteed to be reached given perturbation of the seed set, as long as there are no other interfering signals. This suggests that modules should be somewhat robust to the updating scheme chosen, although additional work is needed to explore module robustness in depth. Pathway module sets, furthermore, provide the logical domain of influence of the seed set [41]. This is a way to estimate which other variables are stabilized (guaranteed to be in a certain state) based on perturbation of a given seed set, regardless of the state of the other variables in the network.
One drawback to this dynamical modularity methodology is that it is still NP-hard in general, and analysis of moderately-sized networks (such as the drosophila parasegment or the T-LGL leukemia network) becomes difficult.
However, it is possible to sample even very large networks with low seed set size and greedy selection criteria to estimate optimal covers. We leave additional analysis of efficient heuristics for finding optimal covers for future work as it is outside the scope of this paper.
Despite the computational difficulty of the problem, our methodology is a step toward understanding how network components interact based on the micro-dynamics of node state updates, and how these dynamical building blocks give rise to the macro-dynamics of network behavior. The authors would also like to thank Deborah Rocha for scientific copy editing.

code availability
The code developed for this paper is made available at https://github.com/ tjparmer/dynamical _ modularity.
[21] Alexander J Gates and Luis M Rocha. Control of complex networks requires both structure and dynamics. Scientific reports, 6:24456, 2016.
[30] Inman Harvey and Terry Bossomaier. Time out of joint: Attractors in asynchronous random boolean networks. In Proceedings of the Fourth         Inputs: The parasegment SPN (size n = 60) extends the single-cell model by including all 14 internal nodes in each of four cells, where WG and HH influence both of the cell's neighbors (periodic boundary conditions assumed), and an additional node, SLP, acts as an input to each cell (see [1] and [2] for a full description).

algorithms
Pathway modules can be discovered by a breadth-first search (BFS) on the DCM (algorithm 1, see Fig. 1), given the caveat of accounting for the temporal dynamics of state updates. Input to the algorithm is a DCM, a seed set S 0 , and a perturbation type (pinning or pulse). The algorithm tracks which s-units and t-units fire in each time step. As with regular BFS units are visited via a first-in first-out (FIFO) queue, but a separate queue is created for each time step. The time counter is initialized at 0 and increases once all of the units in the associated queue have been visited.
T-units are added to the queue associated with the current time step t if their threshold is met, while s-units are added to the queue associated with time step t + 1. If the perturbation type is pinning, any given unit is only visited once (as they are afterwards assumed to fire every time step). If the perturbation type is pulse, by contrast, units that previously fired are allowed to be visited again in future time steps.
The algorithm halts when no new s-units or t-units can be added to a queue (logically, no more units can fire) or the set of s-units and t-units that fire at time t = i is equivalent to the set of s-units and t-units that fired at a previous iteration t = j < i (a cycle is reached) 1 .
The algorithm returns a hash table keyed by time step, where each value is the set of s-units and t-units that fire in the associated time step. The total set of s-units that fire, S, can then be easily derived by finding the superset of s-units from each time step (and ignoring all t-units).
This algorithm can be extended to the general case where each seed s-unit x is associated with a firing schedule, 0 if x is not forced to fire at time t}, rather than only allowing pinning or pulse perturbation. In this case, the external perturbations dictated by the firing schedule take precedence over any units already in the queue. Additionally, the delay associated with units can be generalized; here all s-units take one time step to fire (i.e., they have delay dt = 1); however, s-units can be allowed a higher delay dt > 1 such that it takes several (relative) time steps for a signal to reach an s-unit (for example, this may represent cellular processes that operate on a longer time scale than transcription or translation and could be useful for non-synchronous updating schemes [3,4]).
The set of pathway modules discovered by algorithm 1 can then be reduced to a set of complex modules ( algorithm 2, see Fig. 2). Input to the algorithm is a DCM, a maximum seed set size max_s, a set of s-units seeds, and a perturbation type (pinning or pulse). Algorithm 2 initially finds all pathway modules with seed set size s = 1 whose seed is in the set seeds using algorithm 1 and then removes any subsumed modules to create an initial set of complex modules I 1 . It then iteratively finds the set of complex modules I s>1 with increasing seed set size s = 2, 3, 4...max_s. For each size s, pathway modules are found for all possible seeds sets, where seeds are obtained from the set seeds. This set of potential modules is then reduced to a set of complex modules by removing all modules that are not maximal or that contain non-synergistic submodules. The algorithm halts when the maximum seed set size max_s has been reached and returns the set of complex modules found thus far.
The set seeds can be any set of s-units of interest (for example, the s-units representing the network inputs). If seeds= S, the set of all s-units in the network, then all pathway modules (and therefore all complex modules) will be found; if, however, seeds= Λ 1 , then only complex modules with maximal seeds will be found 2 . This is the maximal seed heuristic which finds all core complex modules in the network.
All analysis and code for this paper, including algorithm 1 and algorithm 2, was written in Python 2.7 with the assistance of the CANA package for finding Boolean network DCMs [5]. The code is publicly available at https://github.com/tjparmer/dynamical _ modularity.  new_seeds={S 0 i ∪ {seed} for i ∈ P, seed ∈ seeds} pathway modules P={} for seed in new_seeds: P ← thresholded_BFS(DCM, seed, perturbation_type) for module p ∈ P: if p is subsumed by m for any module m ∈ P: continue else if µ(S 0 p ) − µ(S 0 ps ) − µ(S 0 p − S 0 ps ) = ∅ for any submodule ps of p: continue else: I ← p s = s + 1 return I Figure 2: Pseudo-code for algorithm 2. Input to the function is a DCM, the maximum seed set size max_s, a set of s-units seeds, and a perturbation type ('pinning' or 'pulse'). At each iteration, new seed sets of size s = k are generated by combining seed sets of current pathway modules of size s = k − 1 with the members of seeds (size s = 1). The pathway modules for these new seed sets are found using algorithm 1; this set of modules is then condensed to a set of complex modules by removing those that are subsumed by other modules or that fail a synergy check. For pinning perturbation, the former condition filters any module p where S p ⊂ S m for any other module m ∈ P, whereas for pulse perturbation the temporal order of when units fire must be checked as well. The function returns the set of complex modules I. In the above pseudo-code, {} represents a set. However, if C and K are used as seeds, then a larger module is seen which includes the extra nodes L and M, |S C,K | = 4 (shown in brown). This synergy is found in two complex modules, |S C,B−0 | = 5 and |S A,K | = 6; however, neither of these are core complex modules because M C and M K are not maximal. Table 1: List of the core complex modules in the drosophila single-cell SPN with pinning perturbation. Modules are grouped by seed set size and ordered by resulting size |S i |. The module labels are colored to correspond to Figure 5 in the main text, as well as supplementary figures 4, 5, and 6; membership in a set of related modules Σ is also indicated in parentheses. The s-units are listed in the temporal order that they are visited within the dynamical unfolding process.        10 , and M 6 , while the latter three groupings have associated modules M 10 , M 6 , and M 9 . The subscript on the s-unit labels denotes the cell of the parasegment where i and j are arbitrary cells chosen such that there is no logical contradiction between variable states and cellular addition uses modular arithmetic. As maximal modules with larger seed sets tend to build on the same set of modules with smaller seed sets, particular seed sets are condensed to variables for readability. Given the periodic boundaries of the parasegment and the similarity of complex modules, a significant degree of redundancy in regards to cellular indices can be seen in the maximal module sets. A seed set of size s 6 is needed to resolve every variable in the network. (a) One of the minimal schemata found to reach an attractor, representing the minimal seeds sufficient to drive the dynamics to that attractor with pinning control. Each row represents a cell in the drosophila parasegment and each column the state of the respective node: black represents ON, white represents OFF, and grey represents a wildcard value (i.e., the node state is redundant for convergence to the attractor under pinning perturbation of the seeds). The color of the row label represents the main module group active in that cell during unfolding (pink represents Σ 1 , red represents Σ 3 , and violet represents Σ 4 , see Table 1). (b) Same as panel a, but for the minimal schema sufficient to reach the wildtype attractor. The minimal configurations found are much smaller than the previous estimates without pinning perturbation found in [2]. . This final configuration is in the basin of attraction of the broad stripes attractor [1], and the network will reach this attractor once the seeds wg and en in cell 1 are no longer pinned.
for input node perturbations. This highlights the usefulness of complex modules in developing a qualitative understanding of control in a system. Perhaps unsurprisingly, the same intracellular modules that lead to attractors in the single-cell SPN lead to attractors in the parasegment SPN, when considering six biologically-relevant attractors from [1]. In particular, M SLP−0,nWG−1 , M nWG−0,nHH−1 , M SLP−1,nHH−1 , and M SLP−1,nHH−0,PT C−1 each play a significant role, appearing in at least four of the attractors (see Table 4).
The logical inferences made during the dynamical unfolding process are also useful in the problem of observability. If a network has reached a fixed point, then each variable state is unchanging (i.e., pinned); thus, observing a set of sensor nodes with corresponding s-units S 0 i guarantees that the module M S 0 i has fully unfolded and all the variable states contained in S i are present in the current state of the network. Observing only a small subset of such states (2-3 in the case of the drosophila single-cell SPN, equivalent to less than 18% of the network size) is sufficient to figure out which attractor the network is in. The size of this minimal observability set is naturally upper-bounded by the size of the minimal driver set. For the drosophila parasegment SPN, the minimal observability sets require only two s-units (3% of the network size) when considering only the six biologically-relevant attractors in Table 4.  Table 5: Comparison of dynamical modularity between genetic regulatory networks. The optimal cover Π * is estimated for each network by maximizing the mean dynamical modularity among covers across seed set sizes, assuming pinning perturbation. Six networks are considered: the example GRN from Fig. 2 in the main text, the drosophila single-cell and parasegment SPNs [1], the Thaliana arabadopsis cell-fate specification network [6], the yeast cell-cycle network [7], and the T-LGL leukemia network [8]. For the first grouping of networks, Π * was calculated considering all possible covers of q 6 core complex modules with seed set size s 5, with higher q values used when computationally feasible. By contrast, Π * for the drosophila parasegment and leukemia networks was estimated using the greedy heuristic mentioned in Section 2.3. For each network the mean dynamical modularity of the optimal cover, D(Π * ), the size of the estimated optimal cover, |Π * |, the average module size within the optimal cover, avg. |S i | ∈ Π * , and the maximum module size within the optimal cover, max |S i | ∈ Π * , is shown. Percentages are calculated based on the number of nodes in the network.  Table 6: Comparison of attractor controllability between genetic regulatory networks. Statistics for the driver sets sufficient to fully resolve each fixed point attractor in the network, assuming pinning perturbation, are shown. Driver sets are equivalent to the seed sets of the minimal pathway modules whose unfolding includes all node states found in the respective attractor. The superset of driver nodes needed (based on the set of s-units sufficient to control to each attractor) is compared to that predicted by feedback vertex set theory [9,10]. See the caption of Fig. 5  Min fixed point seed set size (% of network) 3 (25%) 4 (27%) 3 (18%) 6 (10%) 8 (13%) Table 7: Comparison of core complex modules between genetic regulatory networks. The number of core complex modules, along with the maximum module size, is found per network per seed set size s. Additionally, the minimal seed set size found to fully resolve every variable state (min fixed point seed set size) is shown; note that this fixed point is not necessarily an attractor of the original network, due to the pinning of the seeds involved. It is apparent that variable states in each network are quickly resolved as s is increased when considering the largest complex module. In comparing the larger networks, the leukemia network has more core complex modules than the drosophila parasegment per s value, suggesting more complicated dynamics. Percentages are calculated based on the number of variables (nodes) in the network. See the caption of Fig. 5 for a description of the networks used here.